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Abstract. Stable chaos is a generalization of the chaotic behaviour exhibited by 
cellular automata to continuous- variable systems and it owes its name to an under- 
lying irregular and yet linearly stable dynamics. In this review we discuss analogies 
and differences with the usual deterministic chaos and introduce several tools for 
its characterization. Some examples of transitions from ordered behavior to stable 
chaos are also analyzed to further clarify the underlying dynamical properties. Fi- 
nally, two models are specifically discussed: the diatomic hard-point gas chain and 
a network of globally coupled neurons. 



1 Introduction 

Chaos is associated with an exponential sensitivity of the evolution to tiny 
perturbations in the initial conditions, so that the presence of at least one pos- 
itive Lyapunov exponent is considered as a necessary and sufficient condition 
for the occurrence of irregular dynamics in deterministic dynamical systems 
[1]. In fact, the first observation in coupled-map models of stochastic-like 
behaviour accompanied by a negative maximum Lyapunov exponent came 
as a big surprise [2,3]. In order to highlight the unexpected coexistence of 
local stability and chaotic behaviour, the phenomenon was called stable chaos 
(SC). Although the definition sounds like an oxymoron, in practice, there is 
no logical inconsistency, as the irregular behaviour is a transient phenomenon 
that is restricted to finite-time scales. In spite of this restriction, SC is both 
a well defined and meaningful concept, because the transient duration di- 
verges exponentially with the system size and is therefore infinite in the 
thermodynamic limit. Moreover, the stationarity of SC [3] suggests that it 
can represent an interesting platform for studying non-equilibrium phenom- 
ena. A better understanding of SC can be gained by exploring the analogy 
with the chaotic behaviour exhibited by elementary cellular automata [4], 
another phenomenon that can be formally defined only in the thermody- 
namic limit. In fact, as we clarify in this review, SC is a sort of extension 
of cellular-automaton dynamics to continuous-variable systems. In particu- 
lar, the spreading velocity of localized perturbations, a standard indicator 
used to quantify the degree of chaoticity in cellular automata, proves rather 
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fruitful also to characterize SC. However, in this latter case it is necessary 

to distinguish between finite and infinitesimal perturbations (the latter ones 
cannot even be defined in cellular automata, because of the discreteness of 
the local variable) and it is thereby possible to define two conceptually differ- 
ent propagation velocities. This allows giving a fairly general definition of SC 
as that of a dynamics dominated by "finite amplitude" pertubations [5,6]. 

Altogether, one can express the relevant difference between deterministic 
chaos and SC by referring to the relevant flux of information: while in the 
former context, information fiows from the least towards the most significant 
digits, in the latter, it flows from the boundaries towards the core of the 
system. It would be therefore desirable to develop a general formalism able 
to encompass both phenomena. A promising idea is based on the introduction 
of "non democratic" norms which attribute increasingly small weights to the 
sites that are increasingly "far" from the region of interest. Although this 
approach allows quantifying the spatial information flow, it can be hardly 
extended to account for perturbations that have locally a finite amplitude, 
the analysis of which would require a genuine nonlinear treatment. In fact, 
a tool like finite amplitude Lyapunov exponents (FALEs) [7] appears to be 
more appropriate for characterizing SC, although it is not clear how to go 
beyond the maximal exponent (for the absence of a proper scalar product 
definition in this context). 

As we have mentioned above, in systems with a finite number of degrees 
of freedom, SC is a transient phenomenon. One might therefore think of 
using tools and ideas developed for the characterization of transients such 
as those extensively discussed in the nice review by Tel and Lai [8]. One 
must however distinguish between SC and standard chaotic transients (for 
a seminal paper on the subject, see [9]). In the former case, the maximum 
Lyapunov exponent is positive and formulas such as Kapral-Yorke and Pesin 
relations can be invoked to express some properties of the invariant measure 
in terms of the Lyapunov exponents [10]. In SC, a straightforward application 
of the same formulas yields manifestly useless predictions, as they do not take 
into account the spatial information flow that is the key mechanism of SC. 
Accordingly, one must still heavily rely on direct numerical simulations to 
infer the structure of the invariant measure. Nevertheless, we suspect that 
a possible common property of chaotic transients and SC is the presence of 
a strange repeller. In fact, chaotic transients are almost by definition the 
manifestation of trajectories evolving in the vicinity of a repeller, possibly 
characterized by a small escape rate [8,10]. This property seems to clash with 
the absense of unstable orbits in most of the models exhibiting SC. However, 
such models are also characterized by discontinuities in phase-space and here 
below we argue that their smoothing gives birth to a web of unstable orbits. 
We are thereby lead to conjecture that even though SC is accompanied by 
a negative Lyapunov exponent, its very existence requires the presence of 
topological chaos, i.e. of a finite topological entropy. When and whether the 
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resulting transient dynamics is linearly stable or unstable remains however 
to be clarified. 

The review is organized in the following way. In Sec. 2 we briefly introduce 
the reference models that have mostly used to characterize SC. In section 3 
we properly define SC from the scaling behaviour of the transient length and 
discuss its properties in terms of space-time correlations and fractal dimen- 
sions. Then, in Sec. 4, we discuss the relationship with cellular automata 
by suitably encoding the space-time pattern. In particular, we focus our at- 
tention on the indeterminacy of the next symbol as a way to quantify the 
difference between the original dynamics and that of a suitable deteministic 
automaton rule. We also introduce and estimate the propagation velocity of 
localized perturbations. In Sec. 5, we compare SC with the usual determin- 
istic chaos. This is done by smoothing an otherwise discontinuous coupled- 
map model and studying the dependence of standard indicators such as the 
maximum Lyapunov exponent on the smoothness of the dynamic rule. As a 
result, we identify two thresholds: (i) the first one separates the regions with 
postive and negative Lyapunov exponent; (ii) the second, larger, threshold 
separates the region where finite-amplitude perturbations propagate faster 
than infinitesimal ones, from that where the two velocities coincide with one 
another (which is the signature of a standard chaotic evolution). Moreover, 
we compute the multifractal spectrum of the Lyapunov exponent, showing 
that it has a positive tail even when the average exponent itself is negative. 
In section 6 we discuss various order-to-chaos transitions. In fact, the anal- 
ogy with cellular automata reminds us that such rules are not necessarily 
chaotic. The intrinsic absence of a continuous parameter makes it impossi- 
ble to investigate order-to-chaos transitions in the cellular- automata context. 
As this restriction docs not apply to SC, it makes sense and it is desirable 
to investigate the onset of chaotic dynamics in this latter context. We first 
study a coupled-map lattice, the coupling strength being the relevant control 
parameter. The analysis reveals the existence of a fuzzy transition region, 
where regular and irregular dynamics alternate in a complex manner [11]. A 
simple stochastic model is then introduced to gain some further insight. In 
the new setup, the transition is of directed-percolation type [12]. 

For a long time, SC has been found only in abstract mathematical models, 
characterized by the presence of discontinuities or nearly discontinuous ^ 
evolution rules. This restriction has therefore casted some doubts on the 
physical relevance of this phenomenon. In Sec. 7 we discuss a mechanism that 
can generically lead to discontinuities in physically meaningful contexts. The 
mechanisms requires just the presence of (5-like events such elastic collisions 
between particles or spike emissions by neurons. The "non-commutativity" 
of such events represents a genuine source of discontinuities, which may, in 
turn, give rise to SC. A diatomic hard-point gas and a network of coupled 
neurons are discussed in Sec. 7 as examples of such dynamical systems. The 
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neural network model allows us also to discuss an order-to-chaos transition 

that appears to be of "standard" type (the critical region reduces to a single 
point), although the universality class has not yet been identified [13]. Finally 
the still open problems are briefly summarised in Sec. 8. 



2 Models 

Most of the numerical studies of stable chaos have been carried out in a ID 
lattice of diffusively coupled maps [14], 

Xiit + 1) = (1 - e)f{xi{t)) + |[/(a;i_i(t) + f{xi+i{t))] (1) 

where £ G [0 : 1] is the coupling constant and the map of the interval / is 
piece wise linear. 



PiX + qi < X < Xc 

1- {1 - q2){x - Xc)/f] Xx<x<Xc + ri (2) 
^q2 + P2{x - Xc - v) Xc + r]<x<l, 



where Xc = ( 1 — qi ) /pi . The model [5] is a continuous generalization of the sys- 
tems analysed in [2] and [3], which basically correspond to r/ = 0, i.e. to a two- 
branch maps. The map is continuous because the left and right limits in the 
both connecting points do coincide {f{x~) = f{xf), f{x~ +1]) = +??))• 
Occasionally in this review we speak of "quasi-discontinuous" models, imply- 
ing the presence of large but localized (in phase-space) amplifying regions. In 
this context, this amounts to assuming a small but non-zero width r] for the 
middle branch. In next the section we restrict our analysis to the case 77 = 0. 

Since this model is rather artificial (no specific physical problem lies be- 
hind the choice of /, which has been mostly selected for simplicity reasons and 
for coherence with the seminal paper [15]), we find it convenient to consider 
a second type of model, namely a chain of Duffing oscillators 

Xi = --fXi -x^ + D{xi_i + Xi+i) + {1 + G{t) sm2'Kt/Ti)xi (3) 

where 7 controls the dissipation, D the diffusion between nearby sites and 
G{t) is the modulation amplitude that is periodically switched on and off; 
G = ^ for mod {t, T) < Ti and zero otherwise. As discussed in [4], for T2 
long enough, the Lyapunov exponent is negative, so that the evolution must 
eventually converge towards a periodic orbit, as it indeed does. 



3 Definition and characterization of stable chaos 

Simulations of the above defined map have revealed the existence of long- 
lasting transients followed by a sudden convergence towards some periodic 
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orbit. This suggests, and simulations confirm, that the basin of attraction 
of such orbits is so intricate that the convergence is exponential only for 
distances homogeneously smaller than some threshold 9. However, since there 
exist many different periodic orbits, one cannot estimate the transient length 
by determining the distance from an a priori unknown final state. One can 
nevertheless determine the distance d(t,T) of the configuration {xi{t)} (i = 
1,N) at time t from any previous configuration at time r < i. As soon as 
there exists a r-value such d{tc,T) < 0, we can conclude that d{t + t',T + 
t') will tend to 0, indicating that the dynamics converges towards an orbit 
of period t ~ t. As shown in Fig. 1, the average (over different choices of 
the initial conditions) transient time may increase exponentially with the 
chain length, i.e. with the phase space dimension. This indicates that in the 
thermodynamic limit, the relevant dynamical regime is not the asymptotic 
periodic behavior that is practically unreachable, but what one would naively 
consider a transient regime. This scenario is reminiscent of the disordered 
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Fig. 1. Average transient duration versus the chain length for the diffusively cou- 
pled lattice of maps (2), for pi = 2.7, qi = 0., 77 = 0, g2 = 0.07, p2 = 0.1, and for 
the coupling strength e — 2/3. 

regime in directed percolation, which, in finite systems, has necessarily a finite 
lifetime, as the dynamics sooner or later is absorbed by the homogeneous 
state [16,17]. Irrespective of this difficulty, the disordered regime is a true 
"phase" in the statistical- mechanics sense, as it is stable in infinite systems, 
i.e. when the thermodynamic limit is taken before the infinite-time limit. It 
should, however, be noticed that in SC the "absorbing state" is not just a 
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single homogenous configuration, but may be a set of different and possibly 
exponentially long orbits. 

A second striking character of the transient is that the maximum Lya- 
punov exponent turns out to be negative. Like for the very existence of SC, 
this statement is formally correct only under the assumption of taking first 
the thermodynamic limit. In practice, it is sufficient that the transient du- 
ration is long enough to guarantee a good statistical convergence. From the 
data reported in Fig. 1, one can see that this is not a limitation at all, since 
already in a lattice of 100 maps, the periodic state is practically unreachable. 

The very fact that the transient is Lyapunov-stable makes it substan- 
tially different from the chaotic transients that have been often found and 
attributed to the existence of some chaotic saddle of high dimensionality [8] . 
This is all the way more surprising once we notice that the "transient" dy- 
namics is far from regular. In fact, simulations reveal that both spatial and 
temporal correlations decay exponentially. An example is reported in Fig. 2, 
where we plot 

where (•) denotes an ensemble average 




Fig. 2. Spatial and temporal (smoother curve) correlations for the same parameter 
values of the single map as in Fig. 1 and e = 0.608. 



It is natural to characterize the invariant measure also in terms of its frac- 
tal dimension. Since the whole Lyapunov spectrum is negative, one cannot 
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invoke the Kaplan- Yorke [18] formula to predict the number of active de- 
grees of freedom. Actually, such a formula would imply that the dimension is 
equal to zero and this is in fact true for the asymptotic attractor. Therefore, 
we must rely only on direct numerical computations. More precisely, we have 




Fig. 3. Correlation dimension of the spatial embedding for e — 0.6008 (upper 
panel) and e — 0.608 (lower panel). Dotted, dashed and solid curves correspond to 
embedding dimension e — 1, 2, and 3, respectively. 



decided to compute the correlation dimension [19] of spatial sequences of vari- 
ables [20]. In other words, we have constructed embedding spaces of the type 
Xi(t),Xi+i{t), . . . ,Xi+eit), for e = 1,2,3. In each case, we have counted the 
number of pairs of points Af{e, S) that are separated by a distance larger than 
Zi in a space of dimension e. Afterwards, we have determined the dimension 
as the effective derivative, i.e. 

„ , d\ogJ\f , , 

Formally, the correlation dimension is the limit of D2{e,A) for Z\ ^ 0. As 
for small A, N' is affected by statistical fluctuations due to the finite number 
of points, the relevant question is whether the limiting behaviour sets in for 
distances that are numerically accessible. In Fig. 3 we report the results for 
two different values of the coupling strength, e ~ 0.6008 and e = 0.608, which 
correspond to an ordered and chaotic regime, respectively. Even in the or- 
dered regime, the fractal dimension is finite, as revealed by the plateau, whose 
height is independent of the embedding dimension. The non-zero value of the 
dimension reflects the disordered spatial structure, i.e. the existence of spatial 
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chaos. Therefore, already from this simple case we can conclude on the neces- 
sity to go beyond the standard Lypapunov-exponent analysis. In the chaotic 
regime, the effective dimension is larger and grows with the embedding di- 
mension e (see lower panel in Fig. 3). However, in the absence of theoretical 
arguments, we cannot definitely conclude whether the dimension will satu- 
rate for e 00, indicating the existence of a low-dimensional attractor, or 
whether it diverges, suggesting some form extensivity [20] . 

4 Relationship with cellular automata 

The existence of a stochastic-like dynamics accompanied by an exponential 
contraction of infinitcsimally close trajectories suggests an analogy with the 
so-called chaotic cellular automata (CA) rules [21]. In fact, in a finite lattice, 
any CA rule must eventually produce a periodic orbit, since the number of 
distinct states is finite, namely , where B is the number of states of the 
local variable and N is the number of lattice sites. What makes a chaotic rule 
different from an ordered one is precisely the time needed to cycle through 
previously visited states: such a time is exponentially long in chaotic rules 
[21]. 

A binary representation of the dynamics observed in the the coupled map 

lattice (1) confirms these expectations. The pattern plotted in Fig. 4 are 
indeed very reminiscent of those obtained by iterating CA rules. 





Fig. 4. Two patterns generated by iterating Eq. (1) with the function / defined as 
in Eq. (2) with the same parameter values as in the previous figures and coupling 
strength e = 0.55 (left panel) and 0.7 (right panel). Time flows downwards; black 
corresponds to a;*(t) < 1/2. 



The relationship with CA can be put on more firm grounds, as we dis- 
cuss in the following with reference to the chain of Duffing oscillators (3). 
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The bistable character of the single oscillators suggests a natural way to 

encode the imdcralying dynamics and thereby to explore possible connec- 
tions with CA rules. An appropriate indicator to do so is the indeterminacy 
Ah{'m) of the symbol s' = Si{t + 1), under the assumption that the sequence 
S{m) = {si^r{t), Si-r+i{t), ■ ■ ■ , Si{t), • • • , Si+rCi)}, is obscrvcd at time t (time 
being measured in periods of the forcing term) and where m = 2r + 1. The 
indeterminacy is formally defined as [22] 



where the first sum extends over all sequences of length m generated by 
the chain dynamics, and the second sum to the two values of the sym- 
bol s' . P{S{m)) is the probability to observe anywhere the sequence S{m); 
P{s'\S{m)) is the conditional probability that the observation of the symbol 
s' at time t + I on the site i is preceded by the sequence S at time t in a 
window of length m centered around the site i. When the knowledge of S(m) 
allows to perfectly predict s', then the indeterminacy is zero. In this case, 
the symbolic dynamics is perfectly equivalent to that of a CA defined over 
a window of length m. In Fig. 5 we report the data for T2 = 8, 18 and 20. 
Ah is a non decreasing function of m, since the more we assume to know on 
the past, the smaller must be the uncertainty on the future. If Ah becomes 
exactly equal to for a finite m, then we can conclude that the dynamics 
is perfectly reconstructible from a CA with a finite interaction range. In all 
cases we see that for m larger than 15, the curves saturate revealing the ex- 
istence of a residual uncertainty. This does neither imply that the dynamics 
contains some degree of stochasticity, nor that the model has to include longer 
memory terms. One striking such example was discussed in Rcf. [22], where 
Crutchifield applied this approach to the pattern generated by an elementary 
CA, after it was suitably encoded. The indeterminacy of the encoded pattern 
revealed the presence of a residual uncertainty even though the CA rule is 
deterministic and requires only the memory of one past step, while the en- 
coding is even memoryless. The identification of the "optimal machine" in 
generic cases is a typical example of the hardness of inverse problems. 

The analogy with CA suggests to quantify the degree of chaos also in 
SC by determining the velocity vf of propagation of perturbations. Let us 
consider two configurations that initially differ in the interval [— r, r] and let 
ii{t) {ir{t)) denote the leftmost (rightmost) site where they differ more than 
some threshold. Accordingly, we can define the front velocity as 



Within CAs, a finite spreading velocity is considered as an evidence of chaotic 
behaviour [21]. In fact, this is true also in the context of SC, as it can be seen 
in Fig. 6, where we plot the spreading of an initial difference for parameter 
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Fig. 5. Indeterminacy in tlie cliain of Duffing oscillators (3) with 7 = 0.103, D — 
0.0263, oj = 0.56, Ti = T2/I6 and T2 = 8, 18, 20 (from top to bottom). 

values that correspond to ordered and irregular behaviour. There we see that 
the perturbation spreads only in the latter case (see the right panel). In the 




Fig. 6. Propagation of initially localized perturbations in the coupled map lattice 
for the same values as in Fig. 4. Time flows downwards. 

CA language, the velocity Vf is often named the "Lyapunov exponent" [21]. 
In fact the evolution equation of an elementary CA can be formally written 
as a mapping of TZ^ into itself, 



u\t + l) = F\u\t),u'-{t)) u^'it + 1) = F'-{u\t),u'-{t)) 



(8) 
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where u''{t) = E>o*i(*)2"% = Ei<o*i(*)2* and Si{t) =0,1 (for the 
sake of simplicity wc refer to binary automata). In general and F"^ are 
highly singular functions, but this does not forbids to define a sort of Lya- 
punov exponent from the growth rate of an arbitrarily small perturbation 
S. In practice, if two configurations differ only in the interval [— r, r], we can 
equivalently say that the representation points in TZ^ are separated by a dis- 
tance 6 ~ 2"''. Moreover, if the spatial region where the two configurations 
differ increases with a velocity vj, we can state that the TZ^ distance grows 
as exp[(t!/ log 2)f], thus confirming that, apart from a multiplicative factor, 
the velocity plays the role of a Lyapunov exponent. In coupled map lattices, 
the local variable is continuous rather than binary, but this does not change 
the substance of the argument. 

5 Relationship with deterministic chaos 

The original model where SC has been observed for the first time has a 
peculiarity, namely, the discontinuity of the mapping [3]. As a result, the 
distance betwccin two arbirarily close trajectories can suddenly become of 
order when they find themselves on opposite sides of the discontinuity. 

It is therefore reasonable to study the continuous version of the model, i.e. to 
assume 77 ^ in (2). In the limit 7/^0, the map (2) reduces to the original 
discontinuous system. 

Already at the level of the single map (i.e. without invoking any spatial 
coupling), the introduction of an additional branch may drastically modify 
the structure of the corresponding dynamical system. This is clear in the 
simple case qi = q2 = 0, pi = P2 = 2. For rj — 0, the topological entropy is 
H = In 2, as the map corresponds the Bernoulli shift; however, for any ar- 
bitrarily small, but finite, Tj-value, the appearance of a third branch induces 
a jump to H = In 3. For the parameter values that correspond to the SC 
regime discussed in the previous section, the consequc^iicc of a finite r^- value 
is even more striking, as H is strictly equal to zero for ij — 0, while it is 
finite for = 0+. This can be understood, by performing a slightly noncon- 
ventional symbolic analysis. Let us start by recalling that for the original 
parameter values, there exists a stable period-3 orbit, whose points are or- 
dered as pi < p2 < a;i < X2 < p < 3 < 1. Because of the third contracting 
branch, the intc;rval [,T2, 1] is asymptotically squeezed to a point, so that we 
can identify the leftmost point X2, with the right border of asymptotically 
distinct trajectories. Analogously, the interval [0,pi] is also squeezed to zero 
and we can accordingly interpret pi as the right border of the relevant inter- 
val. Finally [p2,a:i] is also squeezed to zero and can be neglected as well. As 
a result, the relevant dynamics is described by the mapping of /i = [pi,P2] 
and I2 = [xi, ^2]: fill) — I2 and /(/2) = /i U l2- It is easy to show that the 
corresponding topological entropy is the golden mean H = log[(l -|- •\/5)/2]. 
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Therefore, we can at least conclude that the introduction of a finite but arbi- 
trarily small rj induces topological chaos in an otherwise stable environment. 

In the single map, the existence of a fractal chaotic repellor can induce 
long transients only for those trajectories that are carefully selected in the 
vicinity of the repellor itself. It is reasonable to conjecture that in spatially 
extended systems, the stable manifold of the repellor forces generic trajecto- 
ries to follow an intricate arrangement of "channels" before landing on some 
periodic orbits. What are, however, the dynamical properties of the lattice, 
when finite ry- values are assumed? 




Fig. 7. Maximum Lyapunov exponent of tlie map lattice (2) as a function of rj, 
while the other parameter values are the same as in Fig. 1. The results have been 
obtained for N = 200, but are practically independent of the system size. 

First of all, it important to notice that while the topological cnropy jumps 
abruptely to a finite value, the Lyapunov exponent exhibits a smooth be- 
haviour (see Fig. 7), i.e. for sufficiently small r?-values [rj < rj* < 3.10"^) 
it stays negative. This implies that the phenomenon of SC is generic (in 
the mathematical sense), even though the window of existence is (at least 
in this context) rather narrow. More accurate information can be extracted 
by performing a multifractal analysis of the Lyapunov exponent [23,24]. In 
particular, we have computed the probability distribution P(A, t) of the max- 
imum Lyapunov exponent A over a time span t. For sufficiently large t, the 
probability P{X,t) is expected to scale as P{X,t) = exp[— G(A)t] where G 
is a dynamical invariant whose operative definition is obtained by inverting 
this scaling relation, G = — {log P)/t. In Fig. 8 we have plotted the results 
obtained for t = 20 and 40. Even though the most probable and average 
Lyapunov exponent is negative (the spectra refer to 77 = 10~^ < ^*), there is 
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a positive tail, in agreement with the conjectured existence of a web of un- 
stable orbits. The smoothed steps on the right of the maximum correspond 
to the number of times a sample trajectory is actually visiting the expand- 
ing branch. The two spectra do reasonably overlap, suggesting that the time 
t = 40 is already in the scaling regime, although finite-size corrections are 
still large (notice, in fat, that the maximum of G has to be, by definition, 
equal to 0). 



T ' r 




Fig. 8. Multifractal distribution of Lyapunov exponents for rj = 1.10~^, where the 
average Lyapunov exponent is still negative. 

Altogether, SC appears to be somehow complementary to the blow-out 
phenomenon discovered in the study of synchronization transitions [25] . While 
analysing the stability of the synchronization manifold, it has been discov- 
ered that in some circumstances, the corresponding (multifractal) Lyapunov 
spectrum altough mostly confined to the negative semi-axis, may exhibit a 
positive tail. In such a case, one has to go beyond the linear stability analysis, 
because whenever the distance is amplified, nonlinear terms are responsible 
for either bringing the trajectory back towards the manifold or letting it es- 
cape away. In the context of SC, nonlinear terms bring the trajectory back 
towards the "invariant manifold" , although the mechanism is perfectly effi- 
cient only in the infinite dimensional limit. 

In order to clarify the mechanisms by which nonlinearities contribute to 
stabilizing the chaotic dynamics, it is convenient to analyse the propagation 
of perturbations. We start by briefly recalling the concept of convective Lya- 
punov exponents [26] . Given a unidimensional lattice model in the stationary 
regime, let us introduce a 5-like perturbation at time t = in the origin 
i = and imagine to monitor the perturbation amplitude Wi(t). Kaneko and 
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Deissler [26] suggested that 

Wi{t) = exp[yl(i/t)i] (9) 

where A{v) represents the (exponential) growth rate of a perturbation in 
a frame moving with the velocity v. It is a priori obvious that for v = 0, 
one recovers the usual Lyapunov exponent, and that for large velocities one 
has to expect negative growth rates. In fact, A{v) has a typically parabolic 
shape with the maximum in zero. All velocities for which A{v) > corre- 
spond to growing perturbations. The limit velocity for linearly propagating 
perturbations is fixed by the marginal stability criterion A{vl) = 0. Instead 
of determining directly A{v), it is more convenient to exploit the chronotopic 
approach set in [27,28], and formally introduce a perturbation with a spatial 
amplification factor 

Wi{t) = e-f^'uiit) (10) 
In our lattice model, the evolution rule for Ui{t) reads 

Ui{t + l,n) = |e-''7'(a;i-i(t))«i-i(t,M)+ (n) 
(1 - e)f{xi{t))ui{t,,i) + |e''7'(a;i+i(t))wi+i(t,M) 

By iterating this recursive equation with suitable boundary conditions (pe- 
riodic conditions are typically optimal, as they reduce finite-size effects), we 
obtain the chronotopic growth rate A(/x). Altogether, an infinitesimal pertur- 
bation Wiit) = exp[A(/i)i — /ii] with a spatial growth rate /i grows in time 
with an exponent A(^). The evolution of the initially localized perturbation 
is connected to A(^) by a Legendre transform 

A{v)=\{^l)-^l\'(^^) ; v^\' (12) 

In order to determine the velocity corresponding to a given /it-value, it is 
necessary to compute the derivative of A(/i). Since the numerical computation 

of dcrivatcs is always affected by large numerical errors, it is convenient to 
perform a few more analytical steps [27]. By introducing, 

Ui{t, jjL + djji) = Ui{t, i) + Zi{t, iijdn (13) 

in the recursive relation (11), we obtain an equation for the deviation Zi{t, /i), 

Zi{t+1, n ) = |e-'^7'(xi_i(i))(^i_i(i,/x)-Ui_i(t,M)) (14) 

+ (1 - e)f'{xi{t))ui{t,ii) + ^e^'^f'{xi+^{t)){zi+l{t,^x)+Ui+^{t,^x)) . 

The knowledge of Zi and of Ui allows determining A'. In fact, by taking the 
H derivative in the definition of the chronotopic Lyapunov exponent, 

Z t^oo t 
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one obtains 

where • stands for the scalar product. In order to better understand the 
selection process of the propagation velocity, it is convenient to go back to the 
evolution of a single exponential profile Wilt) — exp[A(/i)i — fii]. Its velocity 
is obviously V{fj,) = X/ ii. From the Legendre transform we have that 



dV _1 fdX X\ _ A 
dfj, fi \dii fi J ji 



(17) 



Since the perturbation velocity is identified by the equation yl = 0, we see also 
that it corresponds to the minimum of V{^o). In other words, as long as the 
evolution is controlled by linear mechanisms, the slowest among all possible 
fronts is selected. Let us now turn our attention to fronts delimiting finite 




Fig. 9. Linear {vl, solid curve) and front {vp, dashed curve) velocity versus rj. 
Deterministic chaos exists only for rj > rj* . Beyond rjc, vl — vp. 



perturbations. Since even such fronts must have a leading infinitesimal edge, 
vf will be tii?(/i*) for some fi* . It is hard to imagine that /i* is smaller than 
accordingly, either vp = vl or vp > vl. This scenario is perfectly confirmed 
by the study of the model (2). Solid and dashed curves in Fig. 9 correspond 
to Vp and Vl, respectively. There we see that vp is strictly larger than vl for 
T] < rjc ~ 1.210"'^, while above ry^ the two coincide within numerical accuracy. 
One can also notice that the linear velocity is not defined for rj < if , where 
the system is linearly stable and no propagation of infinitesimal perturbations 
can occur. 
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As discussed in Ref. [6], the mechanism responsible for the finite dif- 
ference between vl and vp is that perturbations of increasing amplitude 
(starting from infinitesimal ones) tend to propagate faster and thereby to 
push the corresponding front. These results are indeed fairly general and not 
just restricted to the model considered in this section [29]. Moreover, this 
phenomenology is conceptually equivalent to that observed in the context 
of front propagation (see, e.g. fironts connecting steady states in reaction- 
diffusion systems [30,31]), that is effectively described by the famous Fisher- 
Kolmogorov-Petrovsky-Piskunov equation [32] . 

Altogether we can conclude by stating that the front velocity proves to 
be a useful indicator to identify the presence of SC (in spatially extended 
systems) from the presence of nonlinear propagation mechanisms that cannot 
be accounted for by linear stability analysis [6]. In such a sense, the results 
in Fig. 9 indicate that SC persists up to the second threshold and not just to 
the first one [5]. 

6 From order to chaos 

Once ascertained that SC is a sort of extension of CA chaos to systems 
characterized by continuous variables, it is natural to investigate the possible 
phase transitions, a question that cannot even be posed in CAs, where all 
variables are discrete. The front propagation velocity vp provides the right 
tool to assess the relative stability of the two phases. Let us, in fact, consider 
two initial conditions: a reference trajectory {x^}, and a perturbed one {xi} 
differing only in a finite interval —L < i < L, where it is randomly set. If 
the interval where v{i) is of order 0{1) increases by eating the region where 
the field was initially equal to zero, we can conclude that the chaotic phase 
is thermodynamically stable. 

In Fig. 10 we plot the results of careful computations performed with the 
coupled map lattice (1) for different values of the coupling strength e. The 
£-range has been selected so as to include both the ordered and the chaotic 
phase. In fact, we see that the front velocity is equal to zero (finite) in the left 
(right) part of the figure. However, these two clearly distinct phases are not 
separated by a point-like transition. We find instead a fuzzy region, where 
chaos and order alternate in a seemingly irregular manner. Is this an evidence 
of the complexity that is sometimes invoked to exist at the edge of chaos? 
Pure numerics alone is not sufficient to provide a convincing answer to such 
a diflacult question. 

An exact formulation and solution of this problem requires to control 
simultaneously two trajectories, a task that is nearly impossible. A simpler 
formulation which can nevertheless help to gain some insight on the transition 
consists in assuming a random evolution for the reference trajectory, and thus 
reducing the problem to that of characterizing the stochastic evolution of the 
difference field Vi{t) [12]. In mathematical terms, this amounts to studying 
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Fig. 10. Front propagation velocity in the coupled map lattice (1,2) for the same 
parameter values as in Fig. 1. 



the equation, 



Vi{t + 1) = (1 - e)wi{t + 1) + ^[wi-i{t + 1) + w^+i{t + 1))] (18) 
where 

/. , IN ivi{t)/v w.p. p = ar] 

W,{t+l)=i if Vi{t)<T] (19) 

\^avi{t) w.p. 1-p 

mit + = "-P- ^ = '^^^^^^ if v.it)>, (20) 

\avi{t) w.p. i- — p 

The stochastic 1 /r] ampUficatioii simulates the effect of visiting the expanding 

interval of the map (2). The amplification saturates to take into account the 
boundedness of the dynamics. This is the only element breaking the linearity 
of Vi{t) dynamics. Moreover, for the sake of simplicity, we assume that a 
uniform contraction rate a (an assumption that is basically equivalent to 
set pi = P2 < !)• At variance with the original deterministic model, here a 
detailed numerical analysis of the parameter space (a, /?), reveals that ordered 
and chaotic phases are separated by a standard phase transition (sec Fig. 11) 
that belongs to the directed percolation (DP) type for small enough values 
of T] [33] and seems to be of multiplicative noise type beyond some critical r] 
value [34]. 
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Fig. 11. Phase diagram of the stochastic model Eq. (20). 



Notice that this stochastic model is even more closely related to the prob- 
lem of synchronization between mutually coupled map lattices (see Refs. [35- 
40] for a more detailed discussion), since the assumption of a stochastic evo- 
lution is appropriate everywhere in parameter space including the critical 
region separating the two phases. 

In the SC context, the DP transition is the most relevant one, as it occurs 
precisely in the regime where the evolution is characterized by a negative 
Lyapunov exponent. DP was introduced and is usually discussed in systems 
where the local variable has just two states: 0, and 1. Moreover, the dynamical 
rule is such that I's cannot spontaneously appear in a sea of O's. This is the 
key difference with respect to the present context, where the variable Vi is 
continuous and thereby the 0-state is never perfectly reached (in finite times). 
It is therefore necessary to introduce a threshold to decide whether the 0-state 
has been reached, with the related problem of having to clarify whether the 
results are truly independent of the threshold. In order to settle this issue, 
wc find it convenient to determine the Finite Amplitude Lyapunov Exponent 
(FALE) [7]. Wc do so by first introducing t{W), the average time needed by 
the field norm 

IKi)ll = ^EKI (21) 

i 

to become for the first time smaller than a preassigned threshold W. 
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The FALE can be thereby defined with reference to a sequence of expo- 
nentially spaced thresholds W„ (W„/W„_i = r < 1) as 



logr 



r(W„+i - Wn) 
In the limit r ^ 1 

-1 



A{W) 



driW) 



dlogW 



(22) 



(23) 



In the further limit ^ 0, yl(M^) reduces to the usual Lyapunov exponent. 
In Fig. 12, we see that the FALE while being almost equal to zero at sufR- 




Fig. 12. Finite amplitude Lyapunov exponent of the stochastic model for two differ- 
ent sets of parameter values, both at criticality: circles refer to A = 0.01, L = 256, 
ttc = 0.6055, squares to A = 0, L = 128, and Oc = 0.6063. 



ciently large scales, becomes equal to the true Lyapunov exponent below a 
certain threshold Wc- Accordingly, since for W < Wc, the dynamics is dom- 
inated by the usual Lyapunov exponent, we can safely conclude that when 
the norm becomes smaller than Wc, the absorbing state will be reached with 
probability one and this solves the problem of an unambigious identification 
of the threshold. Moreover, detailed numerical simulations have revealed that 
Wc decreases faster than 1/L, where L is the system size [33]. In discrete- 
variable systems, the minimal non-zero value that W can meaningfully take 
is 1/L (which corresponds to just one active site). As Wc < 1/L, one can con- 
clude that in this stochastic system the scaling range is even broader than in 
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usual discrete systems. Now a comment about the reason why the linear sta- 
bility analysis may not apply at vanishing distances. In fact, when the above 
defined norm of a vector is small, the field can nevertheless be sporadically 
of 0{1). The behaviour of such bursts may represent an obstruction to the 
validity of the linear stability analysis and this is what tells Fig. 12. 

Finally, we recall that FALE have been employed to characterize single 
maps of the type (2) revealing that for sufficiently small r] and for some fi- 
nite W, the FALE is indeed larger (positive) than the standard Lyapunov 
exponent [29]. Moreover, a generalization of the FALE to a comoving refer- 
ence frame allows to formulate a marginal stability criterion that is able to 
predict the velocity on both cases of linear and nonlinear propagation [29]. 
Moreover, coupled maps (2) with 77 = have been also analyzed by Letz 

6 Kantz [41] who introduced an indicator similar to the FALE (i.e. able to 
quantify the growth rate of non infinitesimal perturbations). This indicator 
turns out to be negative for infinitesimal perturbations and becomes positive 
for finite perturbations. This means that a sufficiently large perturbation can 
propagate along the system due to nonlinear effects. This confirms previous 
observations for marginally stable systems [6]. 

7 More realistic models 

In order to test how general stable chaos is, it is natural to start by asking 
when discontinuities or quasi-discontinuities can be expected to arise in the 
physical world. In fact, we have seen that the source of indeterminacy is the 
sudden amplification of the distance between two nearby trajectories, once 
they fall on opposite sides of a discontinuity. In such a case, no matter how 
small the initial distance is, the separation is suddenly amplified to a value 
of 0(1), that is determined by the size of the discontinuity. 

Before exploring the possible occurence of such phenomena, it is impor- 
tant to stress that the discontinuity we arc referring to is not a discontinuity 
in time of the type associated, e.g., to collisions. A d-like collision induces 
an abrupt change of a variable (specifically, the velocity), but this affects 
the difference between two nearly identical trajectories only for a short (in- 
finitesinal) time lapse, after which the trajectories come close again. This is 
illustrated in Fig. 13, where we have plotted the time evolution of a point-like 
particle bouncing elastically on the floor. In the lower panel the time evo- 
lution of the Euclidean distance is represented: only within the short time 
window between the collisions of the two trajectories with the floor, the rela- 
tive distance becomes of order 0(1). This is at variance with the map lattice 
model (1,2), where the distance, once amplified, remains large. 

In the following two subsections, we illustrate some arguments supporting 
the idea that a natural source of such a type of discontinuities is associated 
with an exchange between non-commuting (5-like events. 




Fig. 13. Separation between two nearly equal trajectories of a point particle collid- 
ing elastically against the floor. The distance 6 = sqrt{Sy^ + Sv'^) is plotted versus 
time. 



7.1 A Hamiltonian model: diatomic hard-point chain 

Before introducing the model, it should be remarked that in Hamiltonian 
systems, conservation of volumes implies that the maximum Lyapunov expo- 
nent cannot be negative; at most, all Lyapunov exponents are exactly equal 
to zero. In fact, the Hamiltonian version of SC is the world of marginally sta- 
ble and yet ergodic models and it often goes under the name of pseudochaos 
(see, e.g. [42]). Here, we are moslty interested in emphasizing the analogies 
with SC and for this reason the diatomic hard point gas (HPG) turns out to 
be rather appropriate also for its relationship with billiard-like models, that 
are often invoked in the analysis of pseudochaos. 

The diatomic HPG is a unidimensional system of point-like particles with 
masses m and M that alternate along a line and undergo elastic collisions. 
In the limit m = M, the model is perfectly integrable, since the velocities of 
the two particles involved in any collision are simply interchanged. Therefore 
there is no mechanism leading to a diffusion in velocity space. However, as 
soon as the masses are assumed to be different, all numerical simulations 
suggest that the dynamics is ergodic. On the other hand, it is easy to convince 
oneself that the maximum Lyapunov exponent is still exactly equal to zero. 
The argument is pretty straightforward [43]: since the collision rule is linear, 

, (m- M)u + 2Mv 
u = 

m + M 

, 2mu — (m — M)v 

V — 

m + M 



(24) 
(25) 
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Fig. 14. Evolution of two nearly equal set of initial condition in a diatomic hard- 
point gas. By slihgly shifting the position of the middle particle, the final velocities 
change abruptly when passing across the three-body collision. 

both real- and tangent-space dynamics follow the same rule. As a result, the 
kinetic energy conservation rriivf = E, where rrii = to or M, depending 
on the parity of i) translates into the conservation of a suitably weighted 
Euclidean norm of the perturbation field, namely, m.iSvf . This means that 
the Euclidean norm of any vector is conserved, irrespective of its direction, 
so that all Lyapunov exponents are equal to zero. 

In the absence of deterministic chaos, which is, therefore, the source of the 
stochastic-like behaviour exhibited by diatomic HPG chains? As illustrated 
in Fig. 14, we argue that the source are the discontinuities occurring around, 
e.g., three-body collisions. Let us consider an initial condition like that in 
the left panel of Fig. 14: it gives rise to a sequence of three collisions, 1 — 2, 
2 — 3, 1 — 2 before the particles separate out. By shifting the position of the 
central particle (this is equivalent to moving the initial Xi variable in the 
CML), we pass to the condition depicted in the right panel, which gives rise 
to the collisions 2 — 3, 1 — 2, 2 — 3. Accordingly, the sequence of two-body 
collisions changes abruptly in correspondence of a three-body collision, when 
the three particles find themselves in the same place at the same time. As a 
consequence of this sudden modification, the three final velocities differ in the 
two cases, as it can be appreciated by comparing the two panels in Fig. 14. 
Only in the limit case of equal velocities, there is no discontinuity, since the 
final set is the same for both sequences. In the former case, when starting 
from the sequence vi, V2, W3, one passes first to V2, vi, v^, then to V2, W3, wi, 
and finally to V3, V2, vi. One can easily verify that the final state is the same 
in the latter case too. 
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Fig. 15. Billiard-like representation of the dynamics of a chain of two particles 
moving in an interval with reflecting boundaries. Horizontal and vertical axes cor- 
respond to the coordinates of the first and second particle respectively. Two slightly 
different trajectories are plotted in each panel until the 5th collision. Left and right 
panels correspond to m = M, and M = 2m, respectively. 

The model dynamics can be further clarified by exploring the analogy 
with billiard models. The connection was first discovered in Ref. [44], where 
the authors considered the slightly different model of hard rods. Here, we il- 
lustrate the relatively simple case of a gas of two particles Pi and P2 bounded 
to move between two fixed barriers Bi and Br, located in x = and a; = 1, 
respectively. The linear position of the two particles can be represented as the 
position of a point-particle in the plane and the constraints < xi < a;2 < 1 
imply that the motion is restricted to the triangular region depicted in Fig. 15. 
Collisions with the two mutually orthogonal triangular edges correspond to 
collisions with either the left or the right barrier, while those with the diag- 
onal correspond to interparticle collisions. Finally, the three angular points 
correspond to the only two possible three-body colHsions, B1P1P2, PiP2Br 
and to the synchronous occurrence of the two-body collisions BiPi a.ndP2Br. 
In the equal mass case, there is a perfect correspondence between hard point 
gas and the triangular billiard. Accordingly, we can invoke the conjecture 
raised in [45] that billiards with rational angles (expressed in tt units) are 
necessarily ergodic. The crucial difference that appears as soon as the two 
masses are assumed to differ from one another is that in the billiard-like 
representation, the mass itself assumes a vectorial character. In particular, 
incoming and outgoing velocities are not mutually symmetric in correspon- 
dence of a collsion with the diagonal. However, the most relevant consequence 
is the appearance of true discontinuities. This is illustrated by comparing two 
nearby trajectories which undergo a different sequence of collisions. In the left 
panel of Fig. 15, which refers to equal masses, we see that the small difference 
in the initial velocity generates a slow linear increase of the mutual distance. 
In the right panel, which refers to M = 2m, the two trajectories, although 
starting from the same initial conditions, drastically separate out and find 
themselves very far apart after as few as 5 collisions (see the arrows) . 
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7.2 Neured networks 

A perhaps more interesting example of a model exhibiting stable chaos is a 
network of leaky-integrate-and-firc ncnrons, where exponentially long tran- 
sients have been identified in various set-ups [46,13,47,48]. By following [13], 
the model dynamics for a network of N neurons can be written as a set of N 
differential equations 

N 

where the connectivity matrix gij is defined as 

if i and j are coupled ^^^^ 
otherwise, 

£i is the number of neurons that are connected to the ith neurons, and G is 
the coupling strength. Here, we will limit to consider the case of inhibitory 
coupling, that, in these notations, corresponds to a positive G value. All 
variables are dimensionlcss and suitably rescaled: the "action potential" Vi <E 
[— oo, 1] whenever reaches the limit value Vj = 1, is reset to and a (5-spike is 
thereby emitted and received by all the connected neurons. The parameter c 
controls the relaxation velocity, while w quantifies indirectly the dependence 
of the effect of the spike on the instantaneous value of the action potential. 

When all connections are active, the dynamics rapidly converges towards 
a stationary state characterized by a sequence of evenly spaced spikes (this 
is a so called splay state [49,50]). In the presence of disorder, such as, e.g., a 
small fraction of randomly broken links, the evolution may signficantly differ, 
depending on the coupling strength G. Below a certain critical value, there 
is still a fast convergence towards an ordered state where the neurons fire 
in a fixed order (in agreement with Jin's theorem [51]); for sufficiently large 
coupling constants, the average (over different realizations of the disorder) 
transient length T is exponentially large with the number of neurons [13]. 
This is illustrated in Fig. 16, where we have plotted the average transient 
for different system sizes: squares and circles correspond to G = 0.5 and 1.8, 
respectively. The solid lines are the result of a linear and an exponential fit, 
respectively. The exponential increase of the transient is a clear indication of 
SC, since at the same time, the maximum Lyapunov exponent (after removing 
the zero exponent corresponding to a shift along the trajectory) is definitely 
negative (as shown in [13]). 

Which is the source of such long transients? In between the spikes, the 
single potentials relax independently towards c (a value that is not reached, 
since c > 1). Therefore, like in the diatomic hard point gas, the evolution 
is piecewise linear and one can derive an analytic expression for the map 
as rigorously done in Ref. [13]. In the absence of jumps between different 
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Fig. 16. Transient length for the neural network model (26) for c = 2, w = 4/7 and 
5% of broken links. Squares (circles) correspond to G = 0.5 and (G = 1.8). The 
solid lines have been obtained by means of a linear and exponential fit, respectively. 

branches, the dynamics would be globally stable; the negativity of the Lya- 
punov exponent is a reminiscence of such a stability. However, like in the 
previous cases, there are discontinuities associated with abrupt changes in 
the firing order of the neurons. Let us indeed consider two neurons i and j 
such that gji = 0, while gij ^ and consider two different initial conditions: 
(i) tii(O) = Wj(0) — e, (ii) Wi(0) = Wj(0) + e. A schematic view of the evolution 
is presented in Fig. 17, where the solid line corresponds to the dynamics of 
the zth-neuron, while dashed and dotted line denote the former and latter 
trajectories, respectively. There, one can see that for times larger than ti the 
two trajectories are separated by a finite distance, as a result of a disconti- 
nuity in the dynamical law. This is due to the dependence of the inhibitory 
effect of a spike on the actual value of v (see the multiplicative factor [u + w) 
in Eq. (26)). Being the size of the discontinuity of the same order of the cou- 
pling strength (©(l/A^)), one might argue that this is negligible for N large 
enough. This is not the case, because it has to be compared with the changes 
induced by the smooth dynamics in between two consecutive spikes that is 
of the same order. Moreover, the distance between the two trajectories is 
even amplified to 0(1) in the time interval [ii, t-2\. As ti — 1\ is, by definition 
of the same order of the interspike interval, this effect too is, in principle, 
nonnegligible. 

More recently, stable and yet irregular behaviour has been reported also in 
the context of a slightly different neural network, where the spike are assumed 
to be received with a finite delay r and the spike effects are independent of 
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Fig. 17. Evolution of two neurons with nearly the same potential and asymmetric 
coupling: neuron i couples to j, but not vice versa. The dotted line corresponds 
to the evolution of neuron i; solid (dashed) curve corresponds to the evolution of 
neuron j in case its action potential is initially smaller (larger) than that of neuron 
j- 

V (see [47]). From the point of view of discontinuities, this latter property 
inhibits a persistent amplification of distances between nearby trajectories. 
Nevertheless, the finite-time amplification mechanism is still present and the 
very fact that long-stable transients have been observed is an indication that 
it lasts enough to yield "avalanches" and thereby to a self-sustained irregu- 
lar behaviour. However, one should also notice that "discontinuities" are a 
necessary but not sufficient condition for the onset of SC. 

8 Conclusions 

In the present Review we have thoroughly discussed the phenomenon of sta- 
ble chaos, a type of irregular behaviour occurring in deterministic systems 
that manifests itself as an exponentially (with the system size) long and sta- 
tionary transient. SC differs from usual chaos in that it is characterized by 
negative Lyapunov exponents, but still reatins some features that are remi- 
niscent of deterministic chaos. In fact, by smoothing out the discontinuities 
present in the most typical SC models, induces the multifractal spectrum of 
the maximum Lyapunov exponent to extend to positive values. This, in turn, 
suggests that topological chaos (i.e. a strictly positive topological entropy) is 
a prerequisite for the observation of SC. However, we have shown that lin- 
ear stability analysis does not to provide a convincing description of relevant 
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properties such as the propagation of finite perturbations. In this respect, a 

promising indicator is represented by the finite amplitude Lyapunov expo- 
nent, although there are conceptual difficulties in extending this approach 
beyond the maximum exponent. 

A further interesting question concerns the transition from ordered be- 
haviour to SC. A detailed numerical analysis of a coupled-map lattice reveals 
the existence of a fuzzy region, where behaviour that is neither strictly or- 
dered nor clearly chaotic has been detected. Is this just a difficulty due to 
strong finite size effects, or this phenomenon hides the presence of a genuinely 
"complex" (uncomputable) evolution? In an almost globally-coupled neural 
network, the transition appears to be a standard point-like phenomenon, 
whose universality class is however still unclear. 

The most important question concerns the generality of SC. All models 
where SC has been observed do possess strong localized nonlinearities that 
may reduce to true discontinuities in phase space. The first models where SC 
has been observed are somehow artificial systems with no direct relationship 
with the physical world. However, the discussion of the diatomic hard point 
gas and of the network of pulse coupled neurons, has contributed to clarify 
that discontinuities may spontaneously emerge in models characterized by the 
presence of non communting "(5-like" events (such as two-body collisions or 
spike emissions). Moreover, since we have seen that SC survives a smoothing 
of the coupled-map model, we may also conjecture that the same holds true 
in these latter contexts, once we assume finite collision times or finite pulse- 
widths. 
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